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The last two decades have witnessed tremendous growth in computational power, the 
development of computational fluid dynamics (CFD) codes which scale well over thousands 
of processors, and the refinement of unstructured grid-generation tools which facilitate 
rapid surface and volume gridding of complex geometries. Thus, engineering calculations 
of 10 7 — 10 8 finite-volume cells have become routine for some types of problems. Although 
the Reynolds Averaged Navier Stokes (RANS) approach to modeling turbulence is still in 
extensive and wide use, increasingly large-eddy simulation (LES) and hybrid RANS-LES 
approaches are being applied to resolve the largest scales of turbulence in many engineering 
problems. However, it has also become evident that LES places different requirements on 
the numerical approaches for both the spatial and temporal discretization of the Navier 
Stokes equations than does RANS. In particular, LES requires high time accuracy and 
minimal intrinsic numerical dispersion and dissipation over a wide spectral range. In this 
paper, the performance of both central-difference and upwind-biased spatial discretizations 
is examined for a one-dimensional acoustic standing wave problem, the Taylor-Green vortex 
problem, and the turbulent channel flow problem. 
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sound speed 

specific total energy, e + Ek 
kinetic energy, U 2 / 2 
specific internal energy 

body force per unit volume vector, [f x , f yi f z ] T 
inviscid flux vectors in the x-, y -, and ^-direction, respectively 
viscous flux vectors in the x-, y-, and z-direction, respectively 
specific total enthalpy, e + Ek + p/ p 

channel half- height, reference length for channel flow problem 
scaled wavenumber, 2 nn/N 

reference length for one-dimensional acoustic and Taylor- Green vortex problems 
Mach number, U/a 
mass flow rate 

subiteration level in Runge-Kutta time-advancement 
number of grid points in a particular direction of a domain 

number of wavelengths in a particular direction of a domain, also time- advancement level 
pressure 

Q-criterion, (R : R — S : S)/2 
heat flux vector, [q xi q yj q z ] T 
gas constant 

rotation-rate tensor, (Vu — Vu T )/2 
strain-rate tensor, (Vu + Vu T )/2 
temperature 
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t time 

U conservative variable state vector 

V velocity magnitude, ||u|| 

u velocity vector, [u,v,w] T 

V volume 

W source vector accounting for body forces 

u,v,w components of velocity vector in the x-, y -, and ^-direction, respectively 

x,y,z cartesian coordinates 

Subscripts 
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Superscripts 

* variable normalized by a characteristic flow timescale 

+ variable normalized by wall friction variables 

I. Introduction 

T urbulence is a critical factor in most aero- and propulsion-related flows. For many years, the industry- 
standard approach to incorporating the effect of turbulence in computational fluid dynamics (CFD) 
calculations has been the Reynolds Averaged Navier Stokes (RANS) modeling approach. RANS implicitly 
time-averages turbulent motion, and models its effect on the mixing of species, momentum and energy 
in a flow field. The approach is economical and often allows a CFD calculation to proceed to steady- 
state. However, it makes several simplifying assumptions, and the results will only be as good as the 
RANS turbulence model. Meanwhile, the last two decades have witnessed rapid and explosive growth in 
computational power, the development of CFD codes which scale well over thousands of processors, and 
the refinement of unstructured grid-generation tools which facilitate rapid surface and volume gridding of 
complex geometries. CFD engineering calculations of 10 7 - 10 8 finite- volume cells have become routine. 
Thus, large eddy simulation (LES), which attempts to directly resolve the unsteady motion of the largest 
scales of turbulence, is increasingly of interest for many fluids engineering problems. 1-3 Additionally, hybrid 
RANS-LES approaches such as Spalart’s Detached Eddy Simulation (DES), 4 which attempt to marry the 
best strengths of both methods, are in active use and development. 

Because it seeks to accurately resolve the unsteady motion of the largest scales of turbulence, LES places 
certain requirements on the numerical method used to discretize and advance the Navier Stokes equations. 
In particular, LES requires high time accuracy and minimal intrinsic numerical dispersion and dissipation 
over a wide range of length scales. With respect to the spatial discretization, central difference methods, 
particularly higher-order ones, can satisfy this requirement. However, the stability of these methods on 
realistic engineering meshes, particular unstructured meshes, can be very problematic. Upwind schemes, 
on the other hand, are typically very robust on engineering meshes and enjoy wide use in many production 
CFD codes. These schemes have been very successful in RANS calculations, but their suitability for LES has 
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ratio of specific heats 
thermal conductivity 
viscosity 

vorticity vector, [lj x , cu y , oo z ] T 
enstrophy, u> • /2 

generic flowfield quantity 
density 

viscous stress tensor 

turbulent kinetic energy dissipation rate 


denotes an initial condition 
denotes a wall friction property 
denotes a flow property at the centerline 

computational indices in the x-, y -, and ^-direction, respectively 
denotes a mean flow property 
denotes a flow property at the wall 
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been arguable for several years. Mittal and Moin 5 compared a 2 nd order central difference and a 5 th order 
upwind scheme on LES of flow behind a circular cylinder at a Reynolds number of 3900, and concluded that 
the upwind scheme significantly dissipated the small scales in the wake. Gamier et al. 6 studied the ability of 
several different upwind shock-capturing schemes to properly resolve near-incompressible and compressible 
decaying isotropic turbulence, and found that while the higher-order shock-capturing schemes did reproduce 
some aspects of turbulence, they were still overly dissipative for this problem. Mossi and Sagaut 7 studied the 
characteristics of several different schemes in DNS and LES of subsonic and supersonic turbulent channel 
flow at a wall- variable Reynolds number of 180. A 2 nd order central difference method did best out of 
the schemes studied, while a 3 th order upwind scheme was again overly dissipative. Park et al. 8 compared 
high-order central and upwind schemes on several different flow problems and reached similar conclusions. A 
common conclusion from these studies is that LES sub-grid scale (SGS) turbulence models, while needed for 
LES using higher-order central difference schemes, are unnecessary for simulations using upwind schemes, 
and in fact typically reduce accuracy at increased cost. 

Although much research in DNS and LES has taken place assuming incompressible fluid flow, the interest 
in this work is in spatial discretizations intended for compressible flow due to the wider range of applicability 
for aerospace flowfields. Ducros et al. 9 developed stable second, fourth and sixth-order skew- symmetric- like 
central difference schemes for compressible flows. The same author, in Ref. 10, developed a flowfield sensor 
suitable for discriminating shocks amid compressible turbulence. Honein and Moin 11 examined the skew- 
symmetric central differencing approach for several different treatments for the conservation of energy and 
entropy. Recently, a large number of additional studies in the literature have investigated the suitability 
of various central difference, upwind-based, and hybrid schemes for LES of compressible turbulence in the 
presence of shock waves. Due to space considerations, only a few will be mentioned here. Larson et al. 12 
found that a hybrid approach using skew- symmetric central difference schemes for the bulk of the flowfield, 
and upwind schemes near shock waves, can be well worth the complexity introduced. Similar conclusions 
were reached in a comprehensive study by Johnsen et al. 13 The recent review by Pirozzoli 14 summarizes 
much of the work to date. 

The work reported in this paper grew out of the author’s work on modeling supersonic film cooling using 
RANS. 15 Subsequent attempts to model the test cases examined in that work using a hybrid RANS-LES 
approach with an unstructured upwind-based production CFD code were problematic, and this study was 
initiated by the author’s attempt to assimilate and understand the above references. The purpose of this 
paper is to evaluate the suitability of several different central difference and upwind inviscid flux schemes 
for a simple one-dimensional acoustic standing wave problem, as well as DNS and LES of the Taylor- Green 
vortex problem and a turbulent channel flow problem. The emphasis in this particular work is on the 
characteristics of the inviscid flux schemes themselves, rather than the coupling of various schemes with SGS 
turbulence models. As such, although the grid resolution in the two turbulence problems ranges from DNS 
levels to fine-grid LES levels, no explicit SGS models are included in the calculations. Additionally, although 
the compressible form of the Navier- Stokes equations are solved, all of test cases occur at low subsonic flow 
speeds, suitable for comparison with well-established results from analysis or the literature. 


II. Governing Equations 

The compressible Navier- Stokes equations for a single chemical species are the governing model for this 
study. The time-dependent differential form of these equations in cartesian coordinates is written as follows: 

dU d(F — F v ) fl(G-G v ) a(H-H v ) 

9t dx dy dz W 

where the state vector, U, and the source vector accounting for body forces, W are given by 
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and the inviscid flux vectors, F, G and H, and the viscous flux vectors F v , G v and H v are given by 
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where the components of the viscous shear stress tensor, r, are given by 
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and the heat flux vector as 
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( 5 ) 


A calorically perfect ideal gas equation of state is used: p = pRT , and T = (7 — 1 )e/R. 7 = 1.4 and 
R = 287 J /Kg-K for all cases studied in this work. The transport properties y and A are calculated using 
the Sutherland formulas for air 16 


y = 1.716 x 10 -5 Pa-s 


T 

273 K 


3/2 273K + 111K 
T + 111 K 


(6) 


A = 2.41 x 1CT 2 W/m-K 


T 

273 K 


3/2 273K + 194K 
T + 194 K 


( 7 ) 


III. Numerical Method 

The Navier-Stokes equations are solved by the finite-difference method in the conventional cartesian 
coordinate system. The inviscid flux terms F, G and H are the portion of the equations that deal with the 
convection of mass, momentum and energy. Historically, the discretization of these terms has been a driving 
issue in CFD research since its inception. Four central difference (denoted as CD-2, CD-4, CD-6 and CD-8) 
schemes and four upwind-biased (denoted as UB-1, UB-3, UB-5, UB-7) schemes are evaluated in this work. 
The number in the name for each scheme denotes the formal order of accuracy. As will be shown below, 
the dispersion (phase) error properties of each central difference scheme is the same as the upwind-biased 
scheme of one less order of accuracy (e.g., CD-4 has the same dispersion error as UB-3). Additionally, an 
upwind scheme based on the Fromm scheme (denoted UF-2) is included, as well as a corresponding central 
difference scheme (CF-2) with the same dispersion error properties. UF-2 should be similar to the linear 
slope reconstructions used in many practical unstructured CFD codes. Finally, a compact-upwind scheme 
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(denoted CU-5) from the work of Pirozzoli 17 is included. Conservative differencing was used such that, at 
node the flux derivatives in the ^-direction is given by 

dx Ax 

where F i+ i - k and F^_i - k are interpolated and formed from the surrounding flowfield points. The same 
process is used for the flux derivatives of G and H in the y- and ^-directions, respectively. 

We now turn to specifying how the flux values at the half-nodes are formed for central differencing. Note 
that the inviscid flux terms F, G and H contain combinations of several different primitive flow variables. In 
divergence form central differencing, the flux values themselves are directly differenced. In skew-symmetric 
form central differencing, the flux is formed from a combination of central differences of the flux values 
and the component flow property variables themselves. In addition to the work in Refs. 9 and 11, this 
approach has been investigated for compressible flow by several researchers in recent years. Kennedy and 
Gruber 18 developed a skew- symmetric central difference scheme intended to be suitable for strong density 
variations. Subareddy and Candler 19 developed a skew- symmetric scheme which maintained consistency 
between the kinetic energy in the momentum and energy fluxes. An especially attractive feature of all these 
schemes is that they preserve kinetic energy, resulting in much greater robustness than divergence form 
central differencing. Pirozzoli in Refs. 20 and 21 showed how several different skew- symmetric schemes could 
be recast in a conservative differencing form. Assuming constant j and k (not shown), and discretization in 
the i-direction, his methodology yields 


CD-2: F i+ i = F avg (i,i + 1) (9) 

CD-4: F i+ i = -F avg (i,i + 1) — - [F avg (i — M + 1) + F avg (M + 2)] (10) 

F^+ 1 = ~F avg (i, i + 1) — — - [F avg (i — 1, i + 1) + F avg (i, i + 2)] 

CD-6: j 10 (11) 

T 7^ [F avg (i 2, i + 1) + F avg (i 1, i + 2) + F av g(i, i T 3)] 

16 4 

F^+ i = Yq ^ av g(b ^ "h 1) — — [F avg (i — 1, i + 1) + F avg (i, i + 2)] 

CD-8: + [F avg (i — 2, i + 1) + F avg (i — 1, i + 2) + F avg (i, i + 3)] (12) 

_ 140 ^ avg ^ — 3, i + 1) + F avg (i — 2, i + 2) + F avg (i — 1, i + 3) + F avg (i, i + 4)] 

CF-2: F i+ i = — F avg (i, i + 1) — — [F avg (i — 1, i + 1) + F avg (i, i + 2)] (13) 


where F avg (il, i2) is formed from a symmetric average of all of the individual components of F. As shown in 
Ref. 20, the form of F avg (il,i2) which reproduces the skew- symmetric scheme of Kennedy and Gruber 18 is 


F avg (il, ^2) — ^ (Pn T Pi2) 2 (^hi T ^ 12 ) 
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(14) 


It is interesting to note that if F avg (il, i2) = ^(F^i+F^) then the classic divergence form central difference 
formulas result, though they can be simplified considerably from the form given in Eqs. (9)-(13). The 
divergence form was tested in the course of this study to confirm the well-known instability characteristics 
of the method. Additionally, other skew-symmetric schemes were tested satisfactorily, including a form 
equivalent to that of Ref. 19. However, unless otherwise noted, all central difference results reported in this 
paper use the Kennedy and Gruber flux in Eq. (14). The flux values for G^^.i k and j fc+ i, formed at 
the half-nodes in the j- and /^-direction respectively, are formed in the same way. 

Upwinding can be implemented by various methods in a CFD code. In addition to the work in Ref. 17, 
high order upwind schemes for compressible flow have been developed by Adams and Shariff, 22 Zhong, 23 
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Gerolymos et al., 24 and Rehman et al., 25 to name only a few. In this work, first a leftward-biased interpolation 
for the primitive flow variables (j) = p, u, v, w, and p at the half- node is formed. Again assuming constant j 
and k (not shown), and discretization in ^-direction, then 


UB-1: 0* + r, L = 0* (15) 

UB-3: cj > iH >L = ( (fri—i + 5 0* + 20* +1 )/6 (16) 

UB-5: 0 i+ 1 L = (2^_ 2 - 130*_i + 470* + 270* +1 - 30* +2 )/6O (17) 

UB-7: 0 i+ i jL = (-60*_ 3 + 5O0*_ 2 - 2O20*_j + 6380* + 4280* + i - 760* +2 + 80* +3 )/84O (18) 

UF-2: 0* + i ?L = (— 0i-i + 40* + 0*+i)/4 (19) 

CU-5: 30 i _i L + 60 i+ i ?L + 0i+|,L = (0i-i + 190* + 1O0 *+i)/3 (20) 


The remaining state variables are calculated using the perfect gas law, and the corresponding conservative 
variables U i+ i and the inviscid flux term F i+ i are then calculated from these components. The process is 
then repeated for a rightward-biased interpolation for the half-node using a reflected stencil. For example, 
0*_i_ i R = (20* + 50*+i — 0*+ 2 )/6 for UB-3, and (p i _i R + 60 i+ i R + 30 i+ g ^ = (100* + 190*+i + 0* +2 )/3 
for CU-5. The final value of the inviscid flux at the half- node is provided by the Roe approximate Riemann 
solver, 26 using these two separate left and right states. 
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( 21 ) 


where A roe is the Roe- averaged flux Jacobian matrix. 

The viscous terms are discretized using 2 nd order central differences. For example, in forming F vi+ i j 
we need to calculate r xx , r xy , r xz and q x at the half- node i+ \,j,k. The x-derivative du/dx at i + \,j,k 
is calculated from (iq + i — Ui)/Ax, and similarly for dv/dx , dw/dx and dT/dx. The cross derivative du/dy 
at i + j, A: is found by averaging the derivative computed at the i and i + 1 nodes, du/dy = — 

+ — ^*+iy-i,/c)/4A?/, and similarly for du/dz , dv/dy and dw/dz. The transport coefficients 

at i + |, j, k are also obtained by a simple average of the values at the i and i - hi nodes. The viscous fluxes 
Gvi j+±,k an d H v y ifc+ i, formed at the half- nodes in the j- and /^-direction respectively, are calculated in 
an analogous fashion. 

All schemes were advanced using a low-storage 4-stage, 2 nd order accurate Runge-Kutta time advance- 
ment algorithm 
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(22) 


(23) 


IV. Fourier Analysis of Schemes 

Before proceeding to the numerical results, it is instructive to review the fourier properties of the various 
numerical schemes for a simple scalar variable 0 discretized into N points over a periodic one-dimensional 
domain, 0 < x < l. Following Ref. 27, we assume a sinusoidal trial solution <j)(x) = ex.p(tkx / Ax) , where 
l = >/~T and Ax = l/N. For n wavelengths in the domain, the scaled wavenumber k is equal to 27rn/N, 
such that it is equal to i r at the spatial Nyquist frequency. The exact solution to the first derivative is 
d(j)/dx = (tk/ Ax) ex.p(ikx / Ax) . For a single variable, 0, if the central difference schemes (Eqs. (9)-(13)) are 
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substituted into Eq. (8), the resulting first derivative of becomes 


cd - 2: 


CD-4: (0;_ 2 

dx 12Ax V 


i - 1 + 1 - 02+ 2) 


CD-6: ~l!x — 60Ax — 450i_i + 450i + i — 90i + 2 + ^+3) 


CD-8: -&jc ~ 840Ax — 32(j)i - 3 + 1680i_2 — 672^_i + 672<^ + i — 1680i + 2 + 32<^ +3 — 3^+4) 


CF-2: ~dx ~ 8 Ax ^ l ~ 2 _ ^ i_1 ^i+i — ^+2) 

Similarly, the upwind schemes (Eqs. (15)-(20)) substituted into Eq. (8) result in 

ub - i: f 

UB-3: ~^“T = (.4*1—2 — + 3</>i + 2(^>j+i) 

UB-5: ~^7 = 60Acc ^~ 2( f >i ~ 3 — 60^>j_i + 20 fa + 30</>j+i — 3</>j+2) 


UB-7: 
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(6<f>i—4 — 56<fri-3 + 2520j_2 — 840</>i_i + 2100* + 504<^_|_i — 84<pi+2 + 8^+3) 
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(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


Using standard fourier analysis, the numerical approximation to the first derivative becomes d<p/dx = 
[i{k' r — bk'j)/ Ax] exp(ifcx/Ax), where k' r and k\ are the real and imaginary parts of the modified scaled 
wavenumber, respectively. The real part k' r is associated with dispersion (phase) error, and for the above 
schemes is given by 


UB-l,CD-2: k' r 

UB-3,CD-4: k' r 

UB-5,CD-6: k' r 
UB-7, CD-8: k' r 
UF-2,CF-2: k' r 

CU-5: k' r 


sin(k) 

4 1 

- sin(fc) — - sin(2/c) 

3 3 1 

2 s[n ( k ) ~ ^ sin(2fc) + — sin(3fc) 

^ sin(fc) - 4 sin(2fc) + 4^ sin(3fc) - 4^ s in(4fc) 

3 1 

- sin (k) — - sin(2 k) 

28 1 \ / 2 \ 1 /18 1 

— sin(fc) + — sin(2A:) J ( 1 + - cos 0 ) ) + 3 sin ( fc ) ( 2 “ 18 C0S ^ ~ ig cos ( 2k ) 

2 \ 2 /l N 2 

1 + - cos(fc) J + ( - sin(fc) 


(35) 

(36) 

(37) 

(38) 

(39) 


(40) 


Note that, as shown by Li, 28 k! r for each of the UB schemes is equal to that of the CD scheme of the next 
higher order of accuracy. The imaginary part of the modified scaled wavenumber, &•, is associated with 
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dissipation error, and is zero for the central difference schemes. For the upwind schemes k[ is given by 

UB-1: k[ = 1 — cos (fc) 

12 1 

UB-3: k\ = - — - cos (fc) + - cos( 2 fc) 

Z o O 

11 2 1 
UB-5: k[ = - - - cos(fc) + — cos(2fc) - — cos(3fc) 

UB-7: fc' = j - 333 cos ( fc ) + Yq cos ( 2k ) ~ 3 ^ cos(3fc) + T cos(4fc) 

UF-2: k\ = ^ — cos(fc) + I cos(2fc) 

( 2 _ 18 C0S ^ “ 18 cos ( 2fc d ( 1 + 3 cos ( fc ) ) _ 3 sin ( fc ) ( Jg sin ( fc ) + Ys sin ( 2fc ) ) 

^1 + | cos(fc)^ + Q sin (k)^j 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 


Plots of the real and imaginary parts of the modified wavenumber are shown in figure 1. In general, the 
dispersion error (the difference between k' r and k) is reduced for the higher-order schemes. Note, however, 
that the UF-2 and CF-2 (Fromm-based) schemes, while having only 2nd-order accuracy, have dispersion 
characteristics comparable to UB-5 and CD-6. The compact upwind CU-5 scheme has superior dispersion 
characteristics to all others. All of the central difference schemes are free from dissipation error (the difference 
between k[ and zero), while all of the upwind schemes experience an increasing dissipation error at higher 
wavenumbers. The higher-order UB schemes have a significantly smaller dissipation error than the lower- 
order ones. CU-5 has a smaller level of dissipation than UB-7 through most of the wavenumber range, but 
crosses over to a higher value at the high end of the wavenumber range. 


V. Results and Discussion 

The five central difference and six upwind inviscid flux schemes were tested on three different problems. 
First, the fourier analysis characteristics just derived were validated using a one-dimensional acoustic standing 
wave problem. Next, the ability of the schemes to accurately resolve vortex breakdown and transition to 
turbulence was tested using the Taylor- Green vortex problem. Finally, the performance of the schemes in 
resolving wall turbulence was tested in a turbulent channel flow calculation. 




(a) Dispersion Characteristics (b) Dissipation Characteristics 

Figure 1. Fourier analysis of spatial discretization schemes studied in this work. 
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V.A. One-Dimensional Acoustic Standing Wave Problem 


As a simple means of evaluating the ability of the various inviscid flux schemes to resolve high-frequency 
acoustics or flowfield variations, they were tested on a one-dimensional acoustic standing wave problem. This 
problem assumes a one-dimensional periodic computational domain in x, with 0 < x < 2irl, and the reference 
length / = 1 m. N = 128 grid points are used. The domain is initialized with air at po = 101325 Pa, To = 
298.15 K, and a small sinusoidal velocity variation as a function of x and specified number of wavelengths, 
n : u(x,t = 0) = Uo cos(nx/l), where Uo = 0.1 m/s. Under the assumption of isentropic, small disturbance 
inviscid flow, the Navier-Stokes equations can be simplified to the linearized equations of gas dynamics, and 
an exact solution derived as 


u(x,t) = -Uo 


COS 


n(x — dot) 
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+ cos 


n(x + a 0 t)\ 

i J 


(nx\ ( na 0 t\ 
= C / 0 CO* ( T j CO* I — I 


p(x,t) -po = -poa 0 U 0 


COS 


n(x — aot) \ 

l 


cos 


n(x + dot) 

l 


TT . fnx\ . f na 0 t\ 
podoUo sin ) sin I — — I 


(47) 

(48) 


where t is time, po is the initial density (1.18413 Kg/m 3 ), and ao is the initial sound speed (346.117 m/s). 

Simulations of this problem were carried out for all eleven inviscid flux schemes for six different numbers 
of initial wavelengths in the domain. Initial numbers of wavelengths of n = 1,2, 4, 8, 16 and 32, corresponding 
to 128, 64, 32, 16, 8 and 4 grid points per wavelength (PPW), respectively, were tested. The corresponding 
scaled wavenumbers are k = 2i rn/N = tt/ 64, tt/32, tt/16, 7t/8, tt/ 4 and tt/2. The simulations for n = 1,2,4 
and 8 were run at a Courant-Friedrich-Lewy number (CFL) number of unity: At = 27rl/(Nao) = 1.41823 x 
10 -4 s. However, due to concern that this timestep was too large to adequately resolve the acoustics for 
the n = 16 and n = 32 cases, they were run at a timestep equal to 1/16 of the characteristic acoustic 
timescale of the problem, 27rZ/(na 0 ), resulting in At = 7.09116 x 1CT 5 s and At = 3.54558 x 10 -5 s (CFL 
numbers of 0.5 and 0.25), respectively. It is worthwhile to point out that, for this simple acoustic problem, 
the divergence form central difference schemes also ran successfully, and produced results identical to the 
baseline skew-symmetric central difference schemes. 

The results for n = 1, 2 and 4 are not shown in this paper simply because there is little difference between 
all the schemes in those cases, with the sole exception of UB-1, which exhibits noticeable dissipation even 
at n = 1. However, greater differences between the schemes become evident at n = 8, 16 and 32. The 
pressure perturbation time histories at the first pressure antinode in the domain are shown for these cases in 
figure 2. The left column shows the results for the central schemes, while the right column shows the results 
for the upwind schemes. The calculations were run through 10 acoustic timescales for each case. However, 
for clarity, only the results for the first two acoustic timescales are shown. 

As expected, there is no dissipation present in all of the central difference schemes, though the dispersion 
error of the CD-2 scheme, noticeable at n = 8, becomes more pronounced at n = 16, and severe by n = 32. 
The dispersion error of the CD-4 scheme is also pronounced at n = 32, and the error is noticeable for the 
other central schemes as well. Interestingly, the CF-2 scheme has a dispersion error in the opposite direction 
from the other schemes at n = 8 and n = 16, consistent with the dispersion characteristics shown in figure la. 

The dissipative characteristics of the upwind schemes are clearly visible in the right column of figure 2. 
The UB-1 scheme is severely dissipative for n = 8, 16 and 32. Dissipation for UF-2 and UB-3 is just noticeable 
at n = 8, more pronounced at n = 16, and severe by n = 32. The dissipation present in UB-5, UB-7 and 
CU-5 only becomes significant at n = 32, with UB-7 less dissipative than UB-5, and CU-5 even better. 
Consistent with figure la, the CU-5 scheme displays very little dispersion error, even at n = 32. 

Eqs. (47) and (48) can be modified to account for the effect of dispersion and dissipation discussed in the 
fourier analysis section, resulting in 


( r) qp > 
— 
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— -^naot 
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(49) 

(50) 


where k is scaled from n, and k' r and k[ are obtained for the numerical scheme under consideration from 
Eqs. (35)-(46). There is excellent agreement between the numerical calculations and these modified equa- 
tions, suggesting that the timestep and time- advancement scheme were sufficiently accurate to resolve both 
the acoustics and the numerical effects from each scheme. 
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(d) ri = 16, 8 PPW, k = 7r / 4 




(f) n = 32, 4 PPW, k = tt/2 


Figure 2. Comparison of pressure time histories at the first pressure antinode in the acoustic standing wave problem. 
One-dimensional grid is 128 points. Note that PPW stands for grid points per wavelength. Left column: central 

schemes. Right column: upwind schemes 


V.B. Taylor-Green Vortex Problem 

The Taylor-Green vortex problem is a benchmark case which, from a smooth initial condition, simulates 
vortex stretching and transition to turbulence. It was used as a test case for the 1 st International Workshop 
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on High-Order CFD Methods, held at the 50 th AIAA Aerospace Sciences Meeting, January 7-8, 2012, in 
Nashville, Tennessee (http://www.public.iastate.edu/ zjw/hiocfd.html). DNS results have been calculated 
by Brachet et al. 29 and, more recently by van Rees et al. 30 The latter results are used as the basis for 
comparison in this paper. 

The solution is computed on a periodic cube domain, —irl < x,y, z < irl. The flowfield is initialized with 


u = f/osili (y) cos (y) cos (y) 


w = 0 


P = P o + 


PoUj 

16 


2x\ 


cos T i + cos V T 


2 y 



(51) 


where here l = 0.01m, po = 7271 Pa, To = 298.15 K and Uo = 34.6115 m/s. The sound speed ao = 
346.115 m/s. The flow is effectively incompressible, Mq = Uo/clq = 0.1. The initial condition for p is 
calculated by the perfect gas law, assuming fixed T = To. The transport properties are fixed at their values 
at Tq throughout the calculation. The resultant Reynolds number is Re = poU^l/po = 1600. Calculations 
for each scheme were performed on three grid sizes: 256 x 256 x 256, 192 x 192 x 192 and 128 x 128 x 128 
grid points, with uniform grid spacing throughout the domain. The characteristic timescale of the problem 
is l/Uo = 2.8892 x 10 -4 s. The calculations were run using a constant timestep of 0.001 times this timescale, 
2.8892 x 10 -7 s, which even on the 256 x 256 x 256 grid is well below an approximate CFL characteristic 
timescale of Ax/(Uo + ao) = 6.4465 x 10 -7 s (CFL « 0.45). The calculations were run to 20 characteristic 
timescales. 

It should be mentioned that while the baseline skew- symmetric central difference schemes ran successfully 
to completion on all cases, simulations using divergence form central difference schemes typically failed 
shortly before a normalized time of t* = tUo/l = 4. Baseline central difference simulations were also run on 
a 64 x 64 x 64 point grid with the viscous terms removed from the equations. These simulations confirmed 
the kinetic energy conservation property of these schemes. 

The time evolution of the Taylor- Green vortex flowfield can be visualized by showing iso-surfaces of the 
Q-criterion 31 at various times (figure 3). Q > 0 indicates regions of a flowfield in which vorticity dominates 
over strain. Here, iso-surfaces of Q = 0.1(L r o//) 2 = 1.19797 x 10 6 s -2 are plotted. This calculation was 
performed using the CD-8 scheme on the 256 x 256 x 256 grid. Time is normalized by the characteristic 
timescale, t* = tUo/l, and the iso-surfaces are colored with the normalized velocity magnitude, U /Uo- It is 
evident from the figure that the initial large-scale vortex structure at t* ip 0 steadily evolves toward smaller 
vortex structures (t* = 4). By t* = 8, the flow has transitioned to turbulence, and proceeds steadily to 
a more isotropic state. The symmetries in the flowfield, discussed by Ref. 29, are evident from figure 3. 
Also note that the reduction in overall kinetic energy as the simulation progresses is visible in the velocity 
magnitude color scale . 

Time histories for the mean kinetic energy, measured mean kinetic energy dissipation rate, and mean 
enstrophy for the central difference and upwind schemes on the 256 x 256 x 256 point grid are compared 
with the spectral DNS (which was calculated on a 512 x 512 x 512 point grid) of Ref. 30 in figure 4. The 
mean kinetic energy for the domain is given by 


Ek : m 


-L ( p U U dV 

PoV J 2 


(52) 


The mean measured kinetic energy dissipation rate is £ = —dE^^/dt. This value was calculated from 
the mean kinetic energy history using standard 2 nd order central differences for all but the first and last time 
history points (which used one-sided differences). 

The mean enstrophy for the domain is given by 


/ “~? dv (53) 

where the vorticity vector w at each point in the domain was calculated using 8 th order central differences 
(i.e. Eq. (27) is used for determining the velocity derivative components of the vorticity vector). It is 
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(a) t* = 0 


(b) t* = 4 


(c) t* = 8 



(d) t* = 12 (e) t* = 16 (f) t* = 20 

Figure 3. Iso-surfaces of Q-criterion Q — 0 .1(IJq / 1)' 2 showing the time evolution of the Taylor-Green vortex problem 
at Re — 1600, computed with the CD-8 scheme on a 256 X 256 X 256 point grid. Iso-surfaces are colored with the 
normalized velocity magnitude, U/Uq. t* is a normalized timescale, tUo/l 




arguable whether this methodology is appropriate for the lower-order flux schemes. It does, however, lead 
to a consistent comparison across all schemes. 

It is evident from figure 4 that the central difference schemes at this grid resolution are in very good 
agreement with the reference DNS results for all three mean flow properties. The CD-2 scheme, which like 
all central schemes is non-dissipative, but has limited range of phase accuracy, does exhibit worse agreement 
with respect to 5 and than the other central schemes near the peak of dissipation at t* =9, though it 
is still reasonably accurate. There are greater differences visible for the upwind schemes. Considering the 
mean kinetic energy history first, the UB-1 scheme, not surprisingly, rapidly dissipates all kinetic energy in 
the flow. The UF-2 and UB-3 schemes follow the reference E^^ m time history reasonably well, although some 
differences from the reference DNS are visible, while UB-5, UB-7 and CU-5 are in very good agreement. 
With respect to the measured dissipation rate, 6, UB-5, UB-7 and CU-5 again agree very well with the 
reference DNS, while greater differences are evident for UF-2 and UB-3. Note that UF-2 and UB-3 have 
a lower dissipation peak that occurs earlier in time (at t* = 8.45). The measured 6 for UB-1 was so large 
in comparison to the others that for clarity it is not shown in figure 4. The upwind schemes exhibit the 
largest differences from the reference DNS in the mean enstrophy history, Q m . Again, as expected, UB-1 
rapidly dissipates all enstrophy in the flow. None of the other upwind schemes match the reference DNS, 
though UB-7 and CU-5 do the best and UF-2 performs the worst here. These results are consistent with the 
dissipative characteristics shown in figure lb. 

It is noteworthy that UB-5, UB-7 and CU-5, which still miss the history near peak dissipation at 
t* = 9, still match the overall history of Ek,m and 6 so well. For an incompressible flow, the theoretical value of 
5 is directly proportional to D m , 5 = 2/iD m /p 0 - Though compressibility effects are present in this calculation, 
their magnitude is small at this grid resolution. Therefore, the good agreement of UB-5, UB-7 and CU-5 
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Figure 4. Plots showing the time evolution of the normalized mean kinetic energy, Ek,m /V^ 2 (upper panel), normalized 
measured mean kinetic energy dissipation rate, el/V^ (middle panel), and normalized mean enstrophy, Ctl 2 /V 2 (bottom 
panel) for the reference spectral DNS on a 512 X 512 X 512 point grid, and the central and upwind schemes schemes 
on 256 X 256 X 256 point grids. Left column: central schemes. Right column: upwind schemes. 


with the dissipation history suggests that, at this grid resolution, the numerical dissipation exhibited by the 
upwind schemes at the small flow length scales is largely making up for the physical dissipation captured by 
the DNS and central schemes at those scales. Note that even UF-2 and UB-3, with much greater divergences 
from the history, follow the reference Ek,m time history reasonably well (though £ not as well). 

The dissipative characteristics of the lower-order upwind schemes are plainly visible in figure 5, where 
iso-contours of normalized vorticity magnitude ||a;||Z/Z7o are shown on a section of the x = 0 plane at t* = 8. 
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(i) UB-5 (j) UB-7 (k) CU-5 (1) UF-2 

Figure 5. Comparison of iso-contours of normalized vorticity magnitude ||ca||Z/t/o on a section of the x = 0 plane for 
all schemes at t* = 8 on a 256 X 256 X 256 point grid. Contour levels are 1, 5, 10, 20 and 30. Note that the DNS was 
performed on a 512 X 512 X 512 point grid. 


All plots, except for the reference DNS, are extracted from the 256 x 256 x 256 grid solutions. The portion 
of the x = 0 plane shown is from — ttI < y — ttI/2 (shown as the x-axis in the plots) and irl/2 < z < i d 
(shown as the y - axis in the plots). The contour levels shown are 1, 5, 10, 20 and 30. UB-1 has completely 
dissipated the vorticity on this plane by t* =8, while UF-2 and UB-3 show significant smoothing of the 
vorticity compared to the DNS reference. The higher-order upwind schemes and the central schemes all 
show generally good agreement with the DNS, though again CD-2 does not resolve the vortex structure as 
sharply as the other central schemes. 

For comparison, time histories for £ T / e?m , measured 5, and for the central difference and upwind 
schemes on the 128 x 128 x 128 point grid are compared with the spectral DNS in figure 6. At this lower 
grid resolution, the central schemes still do reasonably well at matching the DNS, though there is a greater 
deviation in 5 and Q m than was present on the 256 x 256 x 256 point grid. The CD-2 scheme now exhibits 
an earlier peak to these properties (at t* = 8.62) than the other central schemes. Unsurprisingly, the upwind 
schemes UF-2 and UB-3 show very significant differences compared to the reference DNS. From t* = 0-8, 

for these schemes is much lower than the DNS history, while the measured 5 is much higher, suggesting 
that the numerical dissipation in these schemes is too large, and resulting in the faster decay of E k :Tn than 
the DNS. UB-5, UB-7 and CU-5 again do significantly better, and match the Ek : m history reasonably well. 

As was mentioned above, results for all schemes at a grid resolution of 192 x 192 x 192 were also calculated. 
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Normalized Time, t* = t U 0 / 1 Normalized Time, t* = t U 0 / 1 


Figure 6. Plots showing the time evolution of the normalized mean kinetic energy, Ek,m /V q (upper panel), normalized 
measured mean kinetic energy dissipation rate, eI/Vq (middle panel), and normalized mean enstrophy, fll 2 /Vq (bottom 
panel) for the reference spectral DNS on a 512 X 512 X 512 point grid, and the central and upwind schemes schemes 
on 128 X 128 X 128 point grids. Left column: central schemes. Right column: upwind schemes. 


The results are intermediate between the larger and smaller grid resolutions. 

V.C. Turbulent Channel Flow Problem 

Turbulent channel flow has been extensively studied by experiments 32 and DNS studies , 33,34 and the mean 
flow profile and statistics have been very well characterized. Therefore, it is an excellent benchmark case 
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for evaluating various inviscid flux schemes for DNS or LES of wall-bounded shear flows. Here, turbulent 
channel flow of air at po = 25331 Pa, T 0 = 298.15 K and a mean streamwise velocity u m = 44.44 m/s is 
simulated. The flow is essentially incompressible (M m = 0.127). The computational domain is 0 < x < 2 tt h 
in the streamwise (x-) direction, —h < y < h in the vertical (y-) direction, and 0 < z < irh in the spanwise 
(z-) direction. The channel half-width, h, is 0.01 m. The Reynolds number based on u m and the full channel 
width, 2/i, is Re2h = PmUm2h/ /i m = 14300. The DNS calculation in Ref. 34 was run at Reynolds number 
based on wall variables Re r = p w u r h/ p w = 395, which is equivalent to Re^h = 13800, about 4% less than 
this calculation. Note that the wall friction velocity u T = \Jr xy ^ w / p w . 

Periodic boundary conditions are used in both the x- and ^-directions. Adiabatic, no-slip wall boundary 
conditions are imposed on the bottom and top (y-) surfaces. In order to prevent the viscous losses at both 
walls from slowly “wearing down” the flow, the integrated area-sum of these forces is computed at each time 
step, and added in as a body force term, f x , uniformly over the entire computational volume (f y and f z are 
zero) . This term effectively reproduces the effect of the pressure gradient needed to sustain a real turbulent 
channel flow, and maintains a constant mass flow rate (rh = 8.266 x 10 -3 kg/s) in the simulation. 

Simulations for all schemes were run on four different grid sizes: 128 x 129 x 128, 96 x 129 x 96, 64 x 129 x 64 
and 48 x 129 x 48 grid points. The grid spacing is uniform in x and z, and stretched by the hyperbolic 
tangent function in y. The wall-normal grid spacing at the walls is 0.001/i. The 128 x 129 x 128 grid is at a 
near-DNS level of resolution for this flow (the DNS of Ref. 34 was run on a 256 x 257 x 256 point grid). The 
grids with progressively coarser spacing in x and 2 were run in order to assess the sensitivity of the results to 
grid spacing in those directions. The grid size in the ^/-direction was held constant because of the sensitivity 
of the central difference schemes to grid stretching. The maximum grid stretch rate (%+i — yj)/(yj ~ Uj-i) 
with 129 points in the ^/-direction is 1.081, and all of the central schemes were stable at this stretch rate. 
With 97 points in the ^/-direction the maximum stretch rate is 1.118, which is stable for the CD-2 scheme, but 
not for the other central schemes. None of the central schemes were stable with 65 points in the ^/-direction 
(maximum stretch rate of 1.201). One additional simulation was run on a 192 x 193 x 192 point grid using 
the CD-2 scheme (maximum stretch rate of 1.047). 

The simulations for all the central schemes on all grids were run using a fixed timestep of 7.0685 x 10 -8 s, 
which corresponds to an approximate CFL number of 2.45, based on the wall-normal spacing of O.OOlh 
at the walls and a sound speed ao = 346.115 m/s. The UB-3, UB-5 and UB-7 upwind schemes were also 
stable at this timestep. However, the UB-1, CF-2 and CU-5 schemes were generally not stable and most 
of their cases were run at a fixed timestep 25% smaller, 5.6548 x 10 -8 s (the exceptions being UB-1 on the 
128 x 129 x 128, which for stability was run at an even smaller timestep of 4.4178125 x 10 -8 s, and CF-2 on 
the 48 x 129 x 48 grid which was stable at the original timestep). All of these timesteps are more than two 
orders of magnitude smaller than the characteristic timescale of near-wall turbulence, which from Choi and 
Moin 35 is p w /{PwU x ) ~ 1.0 x 10 -5 s for a well-resolved simulation of this problem. Obviously, the explicit 
Runge-Kutta time- advancement scheme, which is limited in stability by the fine near-wall grid spacing, is 
likely not optimally efficient for this problem. 

The flow is initialized using a power law mean velocity profile, along with unit wavelength (n = 1) 
oscillations in y and 2 , and n = 2 oscillations in x. This divergence-free initial condition was borrowed from 
Moin and Kim 36 



(54) 


where here the initial centerline streamwise velocity, u c ? o = 50 m/s and Uo = 5 m/s. 

This initial condition, as well as a fully developed turbulent channel flow, can be visualized by showing 
iso-surfaces of Q = 0.1 (u m /h) 2 = 1.97491 x 10 6 s -2 at t* = 0 and £* = 30 (figure 7). Again, Q > 0 indicates 
regions of a flowfield in which vorticity dominates over strain, t* is a normalized timescale, t* = tu m /27rh y 
representing the number of streamwise flowthrough times. The iso-surfaces are colored with the normalized 
velocity magnitude, U/u m . This calculation was performed using the CD-8 scheme on a 128 x 129 x 128 
point grid. Further details of this case follow below. 
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(a) t* = 0 (b) t* = 30 

Figure 7. Iso-surfaces of Q-criterion Q = 0.1 (um/h) 2 showing turbulent channel flow at Re = 14300, computed with 
the CD-8 scheme on a 128 X 129 X 128 point grid, a) initial condition, b) statistically stationary state. Iso-surfaces are 
colored with the normalized velocity magnitude, U/um. t* is a normalized timescale, i* = tu rn /2izh, representing the 
number of streamwise flowthrough times. 


The central difference cases were started from the initial condition, and it was evident from the skin 
friction history that the flow rapidly transitioned to turbulence within the first few flowthrough times. 
Based on the integrated boundary flux values on the viscous walls, the flowfields were statistically stationary 
by t* = 15. Individual flowfield states were then saved every tenth of a flowthrough time from t* = 15-30 
(151 individual flowfield states were saved for each case). These flowfields were then time-averaged, as well as 
space- averaged in the periodic directions x and z. For comparison, an additional averaging process was also 
performed in which the two halves of the channel {—h<y< 0 and 0 < y < h) were also averaged together. 
Due to the adiabatic walls, the mean temperature and pressure increased by 0.28% over t* = 15-30. 

The upwind scheme cases on the 128 x 129 x 128 and 96 x 129 x 96 point grids were started from 
the corresponding CD-2 flowfield state saved at t* = 10, while the upwind cases on the 64 x 129 x 64 and 
48 x 129 x 48 point grids were started from the corresponding CD-2 states at t* = 12 and t* = 15, respectively. 
The higher-order upwind schemes (UB-5, UB-7 and CU-5) typically reached a statistically stationary state 
within five flowthrough times, while the UF-2 and UB-3 schemes, particularly on the coarser grids, typically 
took longer. The flow laminarized in some of the upwind cases: UB-1 on all grids, UB-3 on the 48 x 129 x 48 
point grid, and UF-2 on the 64 x 129 x 64 and 48 x 129 x 48 point grids. All other upwind cases achieved 
a statistically stationary turbulent state, and the same averaging process as had been applied in the central 
difference cases was performed over 15 flowthrough times. 

Because relatively little change in the averaged flowfield profiles occurred when both halves of the channel 
were combined together into one profile, all subsequent results shown will depict the combined average. Time- 
and space-averaged flowfield values are denoted by an overbar (i.e. u,p w ,jl w and u T = yff xy ^ w / p w , where 
T X y,w = Uw\du/dy\ w ). Time- and space- averaged flowfield variances and covariances are denoted by a prime 
and an overbar (i.e. u' ,v r ,w' and u'v'). 

The wall variable normalized velocity profiles for all schemes on the 128 x 129 x 128 point grid are shown in 
the top panel of figure 8. The velocity profiles are normalized by the characteristic wall friction velocity, u + = 
uju T , and the distance from the wall is normalized by the wall friction lengthscale, y + = | y — y w \u T pw / Uw 
Profiles are also shown for the laminar sublayer, u + = y + and the log layer, u + = 2.441n(?/ + ) + 5.2. Also 
shown is the velocity profile from the Re r = 395 DNS calculation in Ref. 34. In the middle panel of figure 8, 
profiles of the normalized square root of the variance (RMS) for each velocity component, u ,+ = Vti/ 2 /u T , 
v r+ = /u T , and w ,+ = Vw' 2 ju T are shown. Profiles of the normalized turbulent shear stress, — u'v' ju \ 
are shown in the lower panel of the figure. Observe that all five of the central difference schemes (left column 
of the figure) are in generally excellent agreement with the DNS calculation for all of the averaged flowfield 
profiles at this grid resolution. All of them slightly underpredict the peak of fx /+ , but in general they follow 
the DNS profile quite closely. For reference, the grid spacings normalized by the wall friction lengthscale for 
the central difference schemes are: Ar + = Axu r p w /p w = 20, A y+ = 0.42, A y+ = 16 and Az + = 10. 
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Figure 8. Comparison of time- and space-averaged profiles for the inviscid flux schemes on a 128 X 129 X 128 point 
grid. Left column: central schemes. Right column: upwind schemes 


Turning to the upwind schemes (right column of figure 8), both UB-7 and CU-5 also generally agree very 
well with the DNS profiles at this grid resolution. In contrast, UB-5, UB-3 and UF-2 have progressively worse 
agreement with the DNS profiles. The u + profiles rise increasingly above the log law, while u ,+ progressively 
peaks at a higher value further away from the wall. The peaks of v ,+ and w ,+ progressively decrease. The 
displacement of the profile for —u'v'jv? T for UF-2 and UB-3 is also evident. As UB-1 laminarized the flow, 
it is not shown in this figure. 

The same comparison in averaged profiles is shown for all schemes on the 48 x 129 x 48 point grid in 
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Figure 9. Comparison of time- and space-averaged profiles for the inviscid flux schemes on a 48 X 129 X 48 point grid. 
Left column: central schemes. Right column: upwind schemes 


figure 9. All schemes generally do poorer at this coarse grid resolution. CD-6, CD-8 and CF-2 have a u + 
profile which lies consistently below the log law and DNS profile, while CD-2 and CD-4 are overall closer 
to the log law but have an erroneous slope. As an LES SGS model will typically adjust the velocity profile 
upward, use of such a model with the CD-6, CD-8 and CF-2 schemes would likely improve agreement with 
the DNS, while the CD-2 and CD-4 schemes would not likely benefit. The central scheme profiles of fx /+ , 
v ,+ and w r+ diverge considerably from the reference DNS, and particularly so for CD-2. 

Only profiles for UB-5, UB-7 and CU-5 are shown in the right column of figure 9, as all of the other 
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Normalized Spanwise Grid Spacing, p w u x Az/g w Normalized Spanwise Grid Spacing, p w u x Az/p w 

Figure 10. Effect of normalized spanwise grid spacing, A z + = Azii T p w / on the ratio of the skin friction coefficient, 
Cf, to the corresponding Cf predicted by the Dean relation. Left panel: central schemes. Flight panel: upwind schemes 


upwind schemes laminarized the flow. UB-5 at this grid resolution yields results similar to those for UB-3 
and UF-2 on the 128 x 129 x 128 point grid. UB-7 and CU-5 are significantly better, but still not very close 
to the reference DNS profiles. 

Another way to visualize the characteristics of these flux schemes as the grid resolution is decreased is 
shown in figure 10. The ratio of the average skin friction coefficient, Cf = to the corresponding 

value of the skin friction coefficient calculated by the Dean expression, 32 dean, is plotted as a function of 
wall-normalized grid spacing in the spanwise ^-direction, Az+ = A zu T p w /p w . The Dean expression for skin 
friction was developed from a literature survey of experimental data, and is 

Cf, dean = 0.073i?e^°- 25 (55) 

where Re^h is calculated from p m , u m and /! m for each case, but is typically 14260. c/ ? dean for these cases 
is 6.68 x 10 -3 . Also shown in the figure is the result for the reference DNS calculation. Note that the CD-2 
scheme has an additional data point compared to the others, which results from the single simulation carried 
out on the 192 x 193 x 192 point grid. Simulations that resulted in laminarized flow are not plotted. 

It is evident from figure 10 that the central schemes provide a good prediction of the skin friction for 
Az + < 10. At larger values of A z + , the central difference schemes begin to overpredict the skin friction. 
This is due to the lack of grid resolution for resolving the small turbulent lengthscales near the wall, and 
typically an LES SGS model would be employed in these cases to model the contribution of the small scales 
to the flow. While for most of the CD schemes this upward trend continues on the progressively coarser grids, 
note that the CD-2 scheme begins to underpredict the skin friction on the coarsest grids. The relatively 
poor dispersion accuracy of CD-2 on comparison to the other central difference schemes is likely the reason 
for this. 

Considering the upwind schemes, consistent with the results shown previously, UB-7 and CU-5 also 
appear to provide a good prediction of skin friction for Az + < 10, and even provide a reasonable prediction 
(within 10%) for Az + < 13. The predictive accuracy of even these two upwind schemes rapidly falls off at 
larger grid spacings however. Extrapolating the line for UB-5, it is reasonable to expect good predictive 
accuracy for this scheme for Az + « 5. However, it is not within 10% for any of the grid resolutions studied 
in this work. UB-3 and UF-2 are very poor in predictive accuracy. 

VI. Conclusions 

Five different central difference schemes, based on a conservative differencing form of the Kennedy and 
Gruber skew- symmetric scheme, were tested and compared with six different upwind schemes based on 
primitive variable reconstruction and the Roe flux. These eleven different inviscid flux schemes were tested 
on a simple one-dimension acoustic standing wave problem, the Taylor-Green vortex problem, and a turbulent 
channel flow problem. 
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In general, all of the central difference inviscid flux schemes studied proved to be generally very accurate 
and stable, provided that the grid stretching rate was kept below approximately 10%. At near-DNS grid 
resolutions, these schemes yielded numerical results comparable to reference DNS results computed with 
spectral codes. At coarser resolutions, the need for an LES SGS model became apparent for these schemes. 
Consistent with a fourier analysis of the dispersion characteristics of these schemes, there was a noticeable 
increase in the accuracy of the results moving from second to fourth order accuracy. There appear to be 
clear benefits to moving to moving to higher dispersion accuracy, particularly at coarser grid resolutions. 

The seventh order upwind biased scheme, and the fifth order compact upwind scheme, also performed 
very well at near-DNS grid resolutions. The fifth order upwind biased scheme does not do as well, but does 
appear to be suitable for well-resolved DNS. Second and third order upwind schemes having large dissipation 
ranges appear to be poorly suited for DNS or LES. 
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